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A model statistically similar to the flight data was analyzed using the same 
methods (excluding analog matching), and the Newton-Raphson technique was 
found to yield superior estimates. An approximate Cramer-Rao bound was 
compared with the error covariance matrix of the model and was found to provide 
information about the reliability of the individual estimates obtained. 








DETERMINATION OF STABILITY DERIVATIVES FROM FLIGHT DATA 
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This report compares the results obtained by using the modified Newton-Raphson 
technique of minimization with the results obtained from simplified-equations, 
analog-matching, and regression methods for a representative lateral-directional 
flight maneuver. The proposed method is applied to obtain stability and control 



derivatives for vehicles encompassing the geometric and flight extremes of modem 
aircraft. An example of the application of the method to the longitudinal equations of 
motion of an airplane is also presented. 
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weighting matrix for observation vector (R X R) 



weighting matrix for a priori estimate vector (m X m) 
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partition of matrix relating the control vector to the 
observation vector ((R-P) X Q) 







pitching moment divided by the moment of inertia about 
the Y-axis, rad/sec2 



intermediate or incremental time, 



3 > ^ X 




IQ) 


X 

ft 

ra 

i—i 

■8 

•rH 

> 

0 

+-> 

0 

4-3 

CQ 


0 

0 

a 

0 

0 

CQ 

0 

0 

a 


CQ 
0 
i— i 

•8 

•rH 

0 

0 

> 

0 

4-> 

0 

CQ 

<4-H 

o 

CQ 

0 

> 

•r-l 

4-^ 

0 

> 

•H 

0 

*13 


0 

0 

4-> 

15 

0 

0 

0 


X 

ft 


0 

o 


o 


bo 

0 

0 


0 

0 

0 

cq~ 


03 

0 

0 

0 

cr 

CQ 

i 

4-3 

CQ 

0 

0 

»— h 

0 

0 


0 

O 


CQ 

0 

0 

0 

0 

cr 

CQ 


CQ 
0 
0 
» — I 

0 

0 

•rH 

0 

o 



0 

> 

0 

0 

1 


•rH 

0 

•i— i 
0 



bo 


g 

•|-H 

0 

0 

y< 

r ■■■< 
0 

bo 

0 

0 

0 

CQ 

0 

0 

tc 

0 



0 

0 

0 


0 

> 

0 

0 

CQ 

X-axis, 

ft 

•rH 

o 

0 

•rH 

0 

ft 

0 

O 

0 

0 

0 

CQ 

fj 

0 

«s 

0 

rfl 

0 

<2 



o 

0 

0 

0 

r? 


O 

ew 

<W 

O 

*W 

O 

ft 

.a S' 




0 

O 

• rH 


O 

rX 

rM 

•f-4 

rH 

^ X 

Hx 


0 

4-> 

O 

1—1 

v 

0 

0 

o 

0 

4-> 

4-J 

0 

0 

ft 

v/J 

0 

0 

•rH 

ab 

«S r*> 

as 

+-> 

O 

3 

0 

rH 

0 

/X 

a 

0 

0 

CQ 

bJD 0 

bCXS 

0 

•rH 

f . 

0 

ft 

1 

ft 

0 

S 

CQ 

0 

0 

a 

o 

r2 

*bo 

8 

<4-H 

O 

rS 

'bO 

8 

O 

0 

^H 

be 

§ 

•Si 

2a © 
© 

•Si 
•a +5 
■a. a> 

“ a 
© * 

£ 

0 

a 

0 

0 

O 

.a 

0 

> 

4-3 

CQ 

0 

ft 

0 

O 

0 

0 

*— H 

0 


bO 

0 

0 

0 

o 

0 

0 

0 


§ 

’-0 

O 

0 

ft 

0 

0 

0 

1 

■S 

0 


bJD 

0 

0 

0 

O 

0 

0 

0 


g 

ft 

O 

0 

ft 

0 

0 

0 

r§ 

'S 

0 


bO 

0 

"0 

0 

O 

0 

0 

0 


0 

O 

•rH 

4-> 

o 

0 

1— H 

CM 

0 

0 

rH 

0 

s 

1 

1 

4^> 

oq 

g 

o 


0 

0 

•rH 

o 

ft 

«4-H 

0 

O 

O 

0 

O 

a 

0 

0 


0 

0 

bQ 

0 

u 

a 

0 

*0 


bO 

0 

0 

0 

O 

0 

0 

0 


0 

r-*H 

bo 

0 

0 

,0 

O 

•+-> 

ft 


o 

0 

o 


CQ 

0 

0 

0 

4-3 


■8 

-8 

Si 


0 

0 

0 

a a 
0 0 
0 0 
4-3 -+J 

05 


i 


0 
O 

0 Jg 
0 ^ 
0 «H 
> O 

0 >> 
43 £ 


0 

O 

«H 


0 

0 

cr 

0 

0 


§ § 
|p? 

<D i 
®'2 


0 

0 

O 


x 

O CD 

0 

ftft 

Pi 

0 


0 

6 

0 

«S 

CQ 

0 

0 

a 

•rH 

4-> 

CQ 

0 

^H 

0 


0 

0 


tH <M 

O Q 

0 „ , 

0 ° 
0 4 -> 
0 


0 

0 

> 


0 


0 


0 
0 0 
4-3 

0r0 

a - 

ft t 

4-J *rH 

CQ 

0 


CQ 

S 


<S1 


o 

b 


CQ. 


CQ 

ft 


0 

ft 


to 


0 

to 


0 

to 


0 

to 


o 

to 


0 


0 


<M cT 

b 


<Mm- 

b 


ft 

(MO 

b 


NW 

b 


estimate for the variance determined with the 
Cramer-Rao inequality 

variance of the estimates from an experimental model 



auxiliary time variable, sec 


X 

•p4 

u 


a 


0 

o 

% 

•r-4 

cd 

> 

O 

o 

u 

o 

u 

u 

0 

0 

ft 

4-3 

u 

a 


o 

ft 

0 
cd 

ft 

1 

u 

'0 


cd 

O 


be 

0 

■ft 

o 

Td 

cd 

U 


o 

0 

CQ 

Td 

Cd 


0 

*bo 

§ 

ft 

§ 

ft 


>> 

d 

§ 

O 1 

0 

>4 

<4h 


O 

o 


o 

0 

ft 


CQ 

0 

d 

•Cj 


*4-3 


£ 


o 

d 

0 

^3 

cd 

& 


WD 

O 


•4-3 

o 

0 

ft 


CO 

0 

d 

ft 


£ 




o 

13 


0 

•i-4 

n3 

cd 


feo 



w >e s- 3 



0 


ft 


0 

a 


o 



d cr S- 
to « 

« ft CQ. 
cd 
to 


0 

d 

cd 

CQ 

0 

£ 

Id 

ft 

*r-< 

o 

d 


0 

ft 

a 

cd 

CQ 


0 

a 


be 

d 

•r4 

4-> 

d 

0 

CQ 

0 

d 

ft 

0 

d 

X 

0 

Hi 


CQ 

u 

*a 

o 

CQ 

d 

0 

ft 

5 


0 

CQ 

O 

ft 

CQ 

d 

cd 

-b 


cd 


d 

cd 

d 

O’ 

4-3 

cd 

ft 

4-3 

*W 

O 

0 
> 
•i— i 
443 

cd 

> 

•r4 

d 

0 

TJ 

0 

0 


Js 

4-3 

02 

0 

4—3 

o 

d 

0 

T5 


d 

cd 

d 

cr 

cd 


d 

0 

> 

o 

4-3 

o 

Td 


d 

o 

'§> 

d 

o 

d 

43 

4-> 

Td 

0 

CQ 

d 


d 

o 

4-3 

o 

0 

> 

cd 

CQ 

0 

4-3 

cd 

o 

•r-4 

"O 

d 

•rH 

0 

ft 

£ 

0 

a 

e cd 

S3 

o 

ffl 


£ 

O 

h-H 

H 

§ 


ft 

ft 

H 

ft 

Q 

ft 

> 

& 

> 

HH 

ft 

ft 

P 

I 

H 


ft 

g 

co 


CQ 

O 

, a 


o 

ft 

0 

d 


u ™ 

o SS 

^4H ft 5 

"d ^ ^4 

v B o 

ft°3! 

3 d -H 
o 73 ft 

Mi 

0 d H 
ft cd 
^ ft . 


cd 

>a 

CQ 

CQ 

♦rH 

d 

o 

u 

u 

0 



CD 




deals with the lateral-directional modes, because this estimation problem is more 
complex than that for the longitudinal modes. After the proposed technique is fully 
developed and evaluated for the lateral-directional modes, an example is given to show 
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REVIEW OF PREVIOUS METHODS 
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To derive an expression for the values of A 
and B that minimizes equation (2), let 
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used in the minimization because the differential equations are never solved for the 
computed state. 
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The method functions used are those suggested by Shinbrot in reference 4 and have the 
form 



where the maximum of j must be greater 
than or equal to the number of unknown ele- 
ments in C. Then define 





By defining 
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Figure 7. Comparison of X-l 5 time histories Once again it is easily shown that this method 

measured in flight and computed by using has the same row independence of C shown 

coefficients obtained from flight data with the hy both the least-squares and the Shinbrot 
Shinbrot method. method. 
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criteria could be proposed, but the quadratic 
Figure 8. Comparison of X-l 5 time histories measured is the most thoroughly analyzed form and 

in flight and computed by using coefficients obtained has the most desirable mathematical char- 

from flight data with the modified least-squares method, acteristi cs . 







j CS ^ 

5 u g * 

: beg p 

* ^ o 

S> o , 
2 -*-> fl ; 

L«w co g 
r o o 

5 w 0 ) cq 
5 d d 

3 .2 « S ■ 

lt%«- 

i^o u : 

< £J . 

i ’O -e «« 

|®§°. 

< ~ -a s; 

C S 03 G 
) .S > 03 

I be § i i 

a o © 

« ^ to © ' 

) S -rH 03 

oS£ 

r o e-i i 

4 £_f 

j a >> . 

> >,'2 X 3 

j tj .2 5 

03 U g , 

< j£ CO .S ! 

I £ o m 

i 5 a ® t 


o o • • • o 


o © • • • o 


— 1 X 

CO CO 


< . 

bl) 


co 2 

rd CO 


H I O •- 

I • l— 1 

_Lj 

co co 


CQ N 

M ~ 

H t> 









The solution for — — can be expressed as 



where the 1 appears in the i tn row. The corresponding partial derivatives with respect 
to constant bias errors in measuring the state variables can be expressed as follows: 


9(z - y) 
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This method is much more efficient than the gradient method because it attempts to 
predict where the local minimum point is and to step directly to it rather than merely 



stepping in the local downhill direction. However, it is much more complex because of 
the computation of the second gradient matrix. This complexity can be reduced signifi- 
cantly by an appropriate approximation to the second gradient matrix which results in 
a method termed either modified Newton -Raphson or quasi-linearization. It is most 
clearly derived by the method of quasi-linearization. 
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where the superscripts i and j are the indices indicating the time sample, and l 
is the total number of samples. 
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Figure 9. Convergence of fit error for the modified Newton-Raphson method. 
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scalar), multiplied by D 2 . For a weighting of zero the a priori values are ignored, a 
for an infinite weighting the flight data are ignored. Shown in the same figure are two 
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The D 9 matrix for weighting wind-tunnel estimates would ideally be the inverse of 




the error covariance matrix of the test data from which the estimates were obtained, 
inasmuch as this would reflect the relative confidence in each of the a priori values. 



onal of D 2 . (These values are shown in table 1. ) rnis nxea tne relative weignung 

among the individual coefficients. The overall weighting of the wind-tunnel coefficients 
to the flight-data-determined coefficients was changed by varying the values of K. 


'O 2 
JMDfi 


CD 

S m 0 

O CD U 
> Oh 

CQ .5 

h *— ■ 

.2 > 8 


! 0 

i* $ *■ 
S u 2 
1 £ ® "S ‘ 

e pC 

> ^ 0 £ r 

! .S o 1 CQ , 
! cd ^ a) 

| 4-4 *r-H 

[■§si 
i S 2 | : 
. 5 J>S. 

) 03 

j £ S3 S . 
! f-c *•-' <D 

l-S 63 : 

1 . 5 

:£*£ ° 
Xh OQ CD 

, *rH m ! 

; ^ © 3 

in $ © 

s M © 

* xn 0 

) 2 © ffl 
< P ea ca’ i 

) U »> 4-3 

^ cd >h d 

2 I *§ © 
j 8 *i ! 

0-3 3 


_H f-H W 

g cd^J o 

DO fl O ' 0 

■ •as a 

: §'’0 fl 

i K S ° 

i 43 m 


•2 | » w 
i 

o j)5° 
fl 5 o £ 

O ^ n fli 


2 «« 
fl J 


i -' i-*-h CD 

S rj S © <D 

:S.®“ 4J £2 

© GO 2 P 


5 0 O 0 

^ rj 4-^> 

^ /- 1 \ /IN 


.2 8 > 
Cd *rH rH 

£3 > «S 

■m o ° > © 

0 W SU T-I S-J 

O ^ n. ^ ^ 

O ^ ^ 0 5 

•rH CCS 

^ ^ L Q) 


<U 42 M 


tb 3 c« © 

•rH -rH 3 


C'S.S 


■c s 

5® o -a 

> Xi g 
, p 43 


,H « S 'O W5 
to o’© Q, fl 1 j 

- 2 ^ -A 2 

a ® S o 


© £ ~ 3 
73 ^ S $ 

>> a eg a 

Tl © T> 

Is a -g » 

£ 2 a p 
ca bJD^ © 




i’S o M 
I 3 S 

2©^ to 

^ S g ,H 

4-> 0 CD 

Ctf O *{Td 4- 1 
0 0 U CCS 

a 

s ° a c 


©Si 

O a 

task's 

I g 

* S|s 
1 S 31 

^ i? w - 

‘Cj U 2 ^ 
£ ctf +3 £; 


S a .2 
J?N 

. O - = 

"O 

.2 © « 
0^0) 
>a 

® g>«, 


O ^ 
liriTJ C 

• rH 0 

0 <3 *-1 

a^s o 

m I a 


T5 0 ^ 

-£3 0 rn 

-4-1 > CQ 

0 > a 

, S © T3 

; 3© 2 

o 'O rt 

! o. wg 
:«i3S 
, M u 

I I rrH Oh 

S 0 & 

I O L cd 

i ^ 0 

! > g ra 

is s * 
8 

| ® fl O 
t ~£ « 

: 2 8i 

> jS ^ • 

J ffi V 

.3 

! g 2^ § 

4 © CJ M 

! 5j a® 

) cr o ^ 
; (» o ? 


^ 8 «* 
o u u 

0 

^ ^ "O 
Eh 0 0 
^ (0 
0 3 

'd ^ w 

0 >s C3 

S 2 ^ 
0 

S .5 

rt o o 

in 

> ^ o 
© a 3 
a ^.2 

4-> 4-> 

coJ O cd 
tj Cd +» 
Ctl H O 
T3 0 

^ -O 0 

fss. 

^00 

I . ^ 

r-H 

ed ,0 rj 

H H CO 

o ^ T3 
§ w 

0 *3 £ 
pO g ctf 
hi r hi 
H K 

O 0 

•rH 4-1 CD 
£»■§ § 
fl 

rt 0 O 

O 

43 m W 

ll 

3 0 

H ^ W 


Five time histories were generated which differed only in the specific noise 





the value closest to the actual value for every coefficient. Although the variance for the 
Newton-Raphson method is generally larger than that for the least-squares and Shinbrot 
methods, all the actual values are included within the ±lcr-p band about the mean. 
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Because the first term in equation (14) is identical to the term derived for the 
Cramer-Rao bound determined under the somewhat idealized assumptions, the approxi- 
mation of the Cramer-Rao bound is available in the computation of the Newton-Raphson 



minimization. Figure 15 compares the standard deviations obtained by using the 
Cramer-Rao bound from the Newton-Raphson computation with those obtained by analyz- 
ing the statistical model discussed in the preceding section and shown in figure 14. It 
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earlier. However, the D 2 weighting matrix was the same as that used previously, 
and only the overall scalar weighting, K, was changed for each vehicle. The overall 
weighting used was that which doubled the mean square error attained with zero weight- 
ing. As in the X-15 data previously analyzed, the a priori values used were the pre- 
dicted values based upon wind-tunnel tests. 




TABLE 2. -VALUES USED FOR THE D x WEIGHTING MATRICES FOR A LIGHT 
AIRPLANE, THE XB-70 AIRPLANE, AND A LIFTING BODY VEHICLE 
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Figure 1 7 . Comparison of light airplane time histories measured in flight and computed by using coefficients 
obtained from flight data with the modified Newton-Raphson method without a priori weighting. K = 0. 
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Large Supersonic Airplane 

The large, supersonic-cruise XB-70 airplane (ref. 14) was chosen to represent 
this extreme in aircraft. The maneuvers selected consisted of an aileron, 5^, maneu- 
ver followed by a rudder, 6 r , maneuver at a later time in the same flight at essentially 
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Figure 20. Comparison of the XB-70 time histories measured in flight and computed by using a priori values. 
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(b) Second maneuver. 
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Figure 22. Comparison of the XB-70 time histories measured in flight and computed by using coefficients obtained 
from flight data with the modified Newton-Raphson method with an a priori weighting which doubles the unweighted 
fit error. K = 80. 
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(b) Second maneuver. 



the three different methods. It is significant that L p , N^, Y p , and 

are approximately the same for all three methods, but L^, and Ng & are rela- 
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In figure 26 flight data are compared with data computed by using the Newton- 
Raphson method with an a priori weighting of K = 22. The fit is still fairly good, but 
some of the maximum excursions fail to match. 



Figure 26. Comparison of lifting body vehicle time histories measured in flight and computed by using coefficients 
obtained from flight data with the modified Newton-Raphson method with an a priori weighting which doubles the 
unweighted fit error. K = 22. 


Figure 27 is a summary of the coefficient estimates obtained by using the three 
methods. For this vehicle L p , Ljg , Njg, Y p , and are approximately the 

same for all three methods. The coefficients N p and appear to remain near the 

unweighted estimates in spite of the weighting. 
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With x =s Ax + Bu and y=fpix + ‘xT u as the constant coefficient differential 



equation, the Newton-Raphson method was used on the so-called short-period longitu- 
dinal aircraft mode where 



The coefficients in G and H are also in A and B. Because only the short-period 
mode was analyzed, only the elements of the first two rows and columns of the A and 
B matrices and the biases on the state variables q and a were variable. 
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Figure 28. Comparison of XB-70 longitudinal time Figure 29. Comparison of XB-70 longitudinal time 

histories measured in flight and computed by using histories measured in flight and computed by using 

a priori values . K = oo . coefficients obtained from flight data with the 

modified Newton-Raphson method without a priori 
weighting. 



large motion in the maneuver, and linear responses. 
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Raphson method and from a priori values. 
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and the superscript i represents the time of the sample, l is the total number of 
samples, and W is the covariance matrix of the noise, n e . 

For the noise, = z e ~ye- the conditional probability density function can easily 
be shown to be 
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and where c T is the true value of the estimate and the superscript i represents the 

time index of the sample. If, in addition, the estimate of c is assumed to be asymp- 
totically unbiased, the Cramer-Rao bound can be written as 
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